#include <cmath>

#include "DFT.h"
#include "CWT.h"

using namespace ypl;


/*!
 * @brief 一维复数序列的连续小波变换
 * @param aIn       -in-  待变换的复数序列
 * @param aOut      -out- 小波变换的结果（复数序列）
 * @param scale     -in-  尺度因子
 * @param waveType  -in-  小波类型
 */
void CWT::operator () (const Array1D<Complex>& aIn,
                        Array1D<Complex>& aOut,
                        Real scale,
                        WaveletType waveType)
{
    if (aIn.dim() <= 0)
        return;

    switch (waveType)
    {
        case WT_HAAR:
        default:
            cwtHaar(aIn, aOut, scale);
            break;
    }
}

/*!
 * @brief 一维复数序列的连续小波变换
 * @param aIn       -in-  待变换的实数序列
 * @param aOut      -out- 小波变换的结果
 * @param scale     -in-  尺度因子
 * @param waveType  -in-  小波类型
 */
void CWT::operator () (const Array1D<Real>& aIn,
                        Array1D<Complex>& aOut,
                        Real scale,
                        WaveletType waveType)
{
    Integral n = aIn.dim();      ///< 数据点数
    Array1D<Complex> acIn(n);
    for (int i = 0; i < n; ++i) {  ///< 实数数组输入转化为复数数组
        acIn[i].real() = aIn[i];
        acIn[i].imag() = 0;
    }

    this->operator()(acIn, aOut, scale, waveType);
}


/*!
 * @brief 一维复数序列的连续 Haar 小波变换
 * @param aIn       -in-  待变换的实数序列
 * @param aOut      -out- 小波变换的结果
 * @param scale     -in-  尺度因子
 */
void CWT::cwtHaar(const Array1D<Complex>& aIn, Array1D<Complex>& aOut, Real scale)
{
    int n = aIn.dim();                  ///< 带变换点数;
    Real delta = (1.0-0) / (n-1);       ///< 样本步长;

    scale = floor(scale/2) * 2;
    if (scale < 2) scale = 2;           ///< 实际 scale 设为不超过 scale 的最大偶数


    DFT dft;
    Array1D<Complex> aSigFt;
    dft(aIn, aSigFt, DFT::Forward);     ///< 对输入信号作傅里叶变换

    Array1D<Complex> aWaveFt(n);        ///< 小波函数的傅里叶变换;
    Real pi = Constants::pi;
    Real coef = 2 * pi / (n * delta);

    int n2 = n / 2;
    for (int k = 0; k < n; ++k) {
        Real wk = (k <= n2) ? coef*k : coef*(k-n); ///< 频率 omega_k;
        aWaveFt[k] = haarFt(scale * wk);
        aWaveFt[k].imag() = -aWaveFt[k].imag();
    }

    Array1D<Complex> &aProduct = aSigFt;
    for (int k = 0; k < n; ++k) {
        aProduct[k] = aSigFt[k] * aWaveFt[k];
    }

    dft(aProduct, aOut, DFT::Backward);
    Real sqrtScale = sqrt(scale);
    for (int k = 0; k < n; ++k) {
        aOut[k] *= sqrtScale;
    }
}


/*!
 * @brief haar小波函数的连续傅里叶变换函数
 *
 * @param  omega -input- “角频率”，与“线频率” f 的关系为： omega = 2*PI*f
 */
Complex CWT::haarFt(Real omega)
{
    Complex i = Constants::i;
    Complex res;

    if (omega == 0) {
        res = 0.0;
	}
    else {
        Real tmp = sin(0.25 * omega);
        res = i * exp(-0.5*i*omega)* tmp * tmp / omega * 4.0;
	}

    return res;
}
